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We study the thermodynamics of a model of pancake vortices in layered superconductors. The 
model is based on the effective pair potential for the pancake vortices derived from the London 
approximation of a version of the Lawrence-Doniach model which is valid for extreme type-II super- 
conductors. Using the hypernetted-chain (HNC) approximation, we find that there is a temperature 
below which multiple solutions to the IfNC equations exist. By explicitly evaluating the free energy 
for each solution we find that the system undergoes a first-order transition between two vortex liquid 
phases. The low-temperature phase has larger correlations along the field direction than the high- 
temperature phase. We discuss the possible relation of this phase transition to the liquid-to-liquid 
phase transition recently observed in Y-Ba-Cu-O superconductors in high magnetic fields in the 
presence of disorder. 



I. INTRODUCTION 

The properties of flux lines in high temperature super- 
conductors in a magnetic field have been under theoret- 
ical and experimental investigation for many years but 
surprises continue to appear. The first-order transition 
line observed in experiments [Q-Q well below the upper 
critical field Hc2{T) is usually interpreted as a melting 
line of a vortex lattice into a vortex liquid. For Y-Ba- 
Cu-0 (YBCO) superconductors, the first-order nature of 
the transition disappears at critical points in both high 
10,^ and low magnetic fields. Unlike the upper crit- 
ical point, which is usually attributed to the effect of 
disorder, the value of the lower critical field seems to be 
determined by the oxygen content of the sample [Q. In 
Ref. a different interpretation of the first order transi- 
tion was given in terms of a transition between two liquid 
phases with different correlation lengths along the mag- 
netic field direction. The disappearance of the first-order 
transition line at low magnetic fields is explained nat- 
urally as a critical end-point in this picture. For fields 
higher than the upper critical point, the effect of disor- 
der is increasingly important. In a recent experiment on 
YBCO a rich phase diagram at fields above the up- 
per critical field was revealed, which included two liquid 
phases separated by a second-order (continuous) tran- 
sition. Although this transition appears in the regime 
where the effect of disorder is large, it is a reversible ther- 
modynamic phase transition the nature of which is yet to 
be understood theoretically. We believe that the results 
of this paper will be a step towards the understanding of 
this liquid-to-liquid transition. 

In Section II we shall introduce a model for the vor- 
tex system in a layered superconductor in a magnetic 
field perpendicular to the layers. We focus for simplicity 



only on the case where the effect of disorder is negligi- 
ble. The model is based on the London approximation of 
the Lawrence-Doniach model and consists of a system of 
pancake vortices in each layer interacting via an effective 
pair potential obtained by integrating out the degrees of 
freedom other than those of the pancake vortices. We 
find that the present model exhibits a first-order transi- 
tion between two liquid phases with very different cor- 
relations along the field direction. The state below the 
transition is very strongly correlated along the field di- 
rection such that the pancake vortices are likely to be 
stacked on top of each other rather as if in a stiff vortex 
line, while above the transition in the high temperature 
phase the vortex line is dissociated resulting in a phase 
with very short correlation lengths. 

We note that our liquid-to-liquid transition differs 
from another phase transition between two vortex liq- 
uid phases, namely the disentanglement transition, that 
has been studied and debated for some time In the 
present model, unlike below the disentanglement transi- 
tion, there is no long range order along the field direction. 
On the other hand, it might be possible to associate the 
present liquid-to-liquid transition with the one observed 
in the recent experiment |^ on YBCO at fields above the 
upper critical point. Since the effect of disorder has not 
been considered in this model, it is impossible to make 
a quantitative comparison of the present result with the 
experimental results in Ref. Also the orders of the 
phase transition are different from each other. We ob- 
serve, however, that a first order phase transition in the 
absence of disorder becomes a continous transition when 
the disorder exceeds a critical amount 0. This will also 
be the case in our model. It is our contention that the 
basic mechanism that drives the liquid-to-liquid phase 
transition in our model is also responsible for the transi- 
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tion observed in Ref. . 

In this paper the hypernetted chain (HNC) approxima- 
tion is used to calculate thermodynamic quantities and 
correlation functions. The HNC approximation is often 
used for describing correlations in classical liquids 
We set up and solve numerically the HNC integral equa- 
tions corresponding to the pair potential for the pancake 
vortex system. The first order transition emerges when 
below a certain temperature the HNC equations admit 
multiple solutions. By calculating and comparing the free 
energy of these solutions we can determine the thermo- 
dynamically stable branch of solutions. We believe it is 
quite rare in the HNC analysis of classical liquids to find 
multiple solutions and first order transtions in this way. 
Calculations on vortex liquids based on the HNC approx- 
imation have also been done in Ref. , but using a pair 
potential which does not include the Josephson coupling 
between the layers as done here. The authors of Ref. ||] 
used only a phenomenological melting criterion without 
considering the possibility of multiple solutions to their 
HNC equations. 

The paper is organized as follows. In the next sec- 
tion we introduce the model, derive the pair potential 
for the pancake vortices, and set up the HNC equations. 
The results based on the numerical solution to the HNC 
equations are presented and discussed in Section III. 



II. MODEL 

A. Derivation of the pair potential 

Our starting point is the Lawrence-Doniach (LD) 
model [|o| for a layered superconductor in a magnetic 
field along the c-axis perpendicular to the layers. With 
the order parameter in the n-th layer denoted by V'n , the 
LD Hamiltonian is given by 
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where do is the layer thickness, d the layer spacing and 
a, /3, ruab, and rric phenomenological parameters. 

In the London approximation the magnitude {ipnl of 
the order parameter '0„(r) = |V'n| exp(i0„(r)) is held 
fixed and the fluctuations of only the phase 0„ (r) is con- 
sidered. We consider the extreme type-II case where the 
magnetic field B = V x A is constant and uniform. The 
Hamiltonian becomes after neglecting constant terms 
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where J = ?i^|i/'prfo/"^a6, a = (27r/$o)A with $0 — 
hc/2e, and 7 ~ {'mc/inabY^'^ is the anisotropy param- 
eter. In the above equation we approximated the inter- 
layer Josephson couphng term as 2(1 — cos(0„+i —(/)„)) ~ 
{4>n+i ~ 0ri)^- It is convenient to introduce two dimen- 
sionless parameters to specify the system. They are 
r = 27r/3J and ry = (a/jd)'^ where /3 = l/fc^T is the in- 
verse temperature and a = (7r(T)~^/2 13 ^jj^e inter-pancake 
vortex spacing with a being the number density of the 
pancake vortices for each layer that is in turn determined 
by the magnetic field B. In terms of the physical param- 
eters, they are given by 
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where A is the penetration depth parallel to the layers. 
High temperatures thus correspond to small values of F, 
low temperatures to high values of F, and r; is a measure 
of the coupling between layers. 

In order to focus on the vortex degrees of freedom, 
we decompose the phase variable 0„ into the "spin- 
wave part" ipn and the vortex part 0„ for each layer, 
(f'n — + ©n, and integrate out the spin-wave part. To 
do this it is convenient to introduce the Fourier trans- 
form: 

'^"W^'TT^Xl / -^:^'fq(y)(i^vii<ind + ik-r). 



dNi 



In this paper we consider a system which consists of a 
stack of Ni layers, on which we impose periodic bound- 
ary conditions in the direction perpendicular to the lay- 
ers; ipn+Ni{ic) ~ ^n(y^)- where q — 2'n:j/Nid and we use 
Ni integer values of j in the sum such that the wavevector 
q G [— tt/c?, Tr/d) belongs to the first Brillouin zone. 0q(k) 
is defined in the same way. For future use we introduce 
the hatted quantity ^„ = (a/d)ipn for each function. It is 
straightforward to integrate over tpq to obtain the vortex 
part of the Hamitonian given by 

J 



n 

^(l-cosgd) 



dNi 



(4) 



(27r)2 fc2 + _cosgd) 



fc2e_,(-k)e,(k) 



The above Hamiltonian represents pancake vortices on 
each layer interacting through some pair potential w„(r) 
which can be read off from the above equation as follows. 
The key point in deriving the pair potential for the pan- 
cake vortices is to note that one can write the pancake 



vortex density pn (r) 
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where 6^(r) is given by d^Qn = -dyQ'„, dyQn = dxQ'n- 

(i) 

In the above expression r„ are the positions of the pan- 
cake vortices. We also note that (V9„)^ — (V6JJ^ and 
—k^Q'q(k) = ^^(k). Finally to get the pair potential 
w„(r) , we write the Hamiltonian in terms of p„ as fol- 
lows: 



For the case where there is no inter-layer coupling, 77 = 0, 
we can see using the above procedure that the pair po- 
tential behaves as 1/fc^, i.e. logarithmically in real space 
and that Tiy describes Ni independent copies of the 
two-dimensional one-component plasma (OCP) where 
charged particles interact with each other via a potential 
which is the logarithm of their spatial separation, all in 
the presence of a neutralizing background charge. The 
long-ranged nature of the interaction is a direct conse- 
quence of neglecting the screening of the magnetic field. 
The neutralizing background is provided here by the ap- 
plied field represented by a. In the presence of the inter- 
layer coupling we can generalize this analysis and we find 
that the pair potential, Wq{'k) = f3vq{k) is given by 
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where — 2TraT. In real space, the potential energy be- 
tween two pancakes in the same layer is logarithmic and 
repulsive: between pancakes in different layers there is 
an attractive Coulomb-like interaction. It is this attrac- 
tive interaction which is the origin of the liquid-liquid 
transition and a low-temperature phase where a pancake 
tends to stack on top of one in the layer below to form a 
stiff vortex line in order to take advantage of the forces 
of attraction between them. 



B. HNC equations 

We shall set up and solve the HNC equations associ- 
ated with the pair potential (^). First the direct cor- 
relation function c„(r) and the pair correlation function 
/i„(r) are related through the Ornstein-Zernike relation. 
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The pair distribution function 5,1 (r) = hn{r) -I- 1 is given 
in terms of the other correlation functions via the HNC 
closure equation M , 



gn{r) = exp [ - /3t;„(r) - c„(r) -I- ft.„(r)] . 
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In the following we shall solve these equations nu- 
merically. It is well known from the previous studies 



pT[ of the two-dimensional OCP using the HNC ap- 
proximation, one should take special care of the long- 
ranged logarithmic potential. In the present model the 
long-range tails in the pair potential are given by the 
terms singular as l/fc^ for small k in {y<j(k). These tails 
behave as (2 — 6qfi)kf~i/k'^. We separate the pair po- 
tential into the short-range part with an arbitrary mo- 
mentum cutoff Q and the remaining long-range part: 



,(k) = 4^'(k) 



(k), where 
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According to Eq. (0) , the long-range part of the pair po- 
tential must be cancelled by a long-range part of the di- 
rect correlation function, i.e. c„(r) — Cn\r) + Cn\r) 



(r) such that Eq. (|^) only involves 



where c„ (r) 
short-range functions: 

5„(r) = exp [ - wi:\r) - c(f)(r) + h,,{r)] . (10) 

Eqs. and ( p^ form a closed set which we solve it- 
erativcly. There are many equivalent iteration methods. 
Here we start from some function Cn\v) and calculate 



= h„(k) - ^a \'k) from Eq. (|). Then we ob- 



tain the updated Cn^ (r) from Eq. ( |10| ) and the definition 
Nn\r). In the numerical calculation it is crucial to have 
rapid convergence. Here we employ the method intro- 
duced in Ref. p^ , where a combination of several pre- 
vious steps were used as an updated Cn \r). We utilize 
the results from the five previous iterations to produce a 
new input function. For small values of F it is sufficient 
to start from Cn'' (r) = as the initial function. We then 
solve the HNC equations for fixed 77 and for increasing 
values of F using the solution obtained for lower values 
of F. 



III. RESULTS AND DISCUSSION 

In this section we present the numerical solutions to 
the HNC equations described above. Our calculation is 
done mainly on a system consisting of Ni — 11 layers. In 
the last part of the section we consider how the solutions 
behave for larger values of Ni . 

One of the most important results of the present work 
is that the HNC equations describing a pancake vortex 
liquid admit multiple solutions over a range of values 
of 77 and F. This arises in the numerical calculation in 
the following manner. We find that the solutions which 
are obtained in the way described above cease to exist 
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above some value Tj. We call these solutions type-I. 
For r smaller than Tj we were able to find another type 
of solutions (type II) which are different from the type- 
I solutions but which approaches the type-I solution as 
r — + F/. The type-II solutions continue to exist as long 
as F > Tjj for some F//. Below F// again the type-II 
solution disappears. Wc then find yet another solution 
(named type-Ill) for F > F// which arc different from 
the type-II and type-I solutions. 

Among the three solutions that exist for F// < F < F/, 
only one will be a true thermodynamically stable solu- 
tion that describes the system in equilibrium. The solu- 
tion which is stable can be determined by finding which 
of the three solutions has the lowest free energy for given 
F. The free energy can be evaluated using the correlation 
functions calculated within the HNC approximation . 
Explicitly we calculate the excess (over the ideal gas) free 
energy F°^, or the dimensionless free energy per parti- 
cle f'' = F°''/2t:JN with N being the total number of 
pancakes in the system as follows: llsl 



/"'^^^E/ rf'r[-c(f)(r)-«;(^)(r) 
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Figure 1 shows a typical plot of the free energy as a func- 
tion of F for fixed -q. We see that the type-II solutions 
are everywhere unstable and that F/ and F// are two 
spinodal points. We can identify Fc where the first order 
phase transition occurs such that above and below this 
value the system is described by the type-Ill and type-I 
solutions respectively. 

It is important to note that this first-order transition is 
not a melting transition nor does it involve the existence 
of long-range order in the direction perpendicular to the 
layers. This can be most easily seen from the structure 
factor defined by 



5'„(k) = 1 + ^ exp(igW)/iq(k). 



Above and below the transition temperature Fc the in- 
plane (i.e. n = 0) correlation is almost the same as can 
be seen from the structure factor in fc-space in Fig. 2 (a) 
and also from the pair correlation function ho{r) in real 
space Fig. 3 (a). We also note that no crystalline long 
range order exists below the transition, since the peaks 
have finite widths. The most dramatic changes are ob- 
served in the inter-layer correlation. This can be seen by 
comparing S'„(k) for n ^ Q above and below the transi- 
tion as in Fig. 2 (b)-(d). Compared to the situation above 
the transition, the system below the transition remains 
highly correlated even when the layers are far apart. As 



one lowers the temperature (increases F) across the tran- 
sition, the pancake vortices that were more or less in- 
dependent of those in neighboring layers become instead 
highly correlated. 




FIG. 1. The excess free energe per pancake f^^ as a func- 
tion of r for 77 = 0.5. The free energy for the type-I (circles) 
and type-Ill (squares) solutions are shown together with the 
unstable type-II solutions (triangles, dotted line). The arrow 
in the middle indicates the transition temperature Fc = 24.5. 
The arrows in the left and right of Fc represent the spinodal 
points Tii ~ 13.0 and Fj ~ 49.0, respectively. 
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FIG. 2. The structure factor S„[K) for (a) n = 0, (b) 
n = 1, (c) n = 2 and (d) n = 5 at 77 = 0.5, where K = ka \s a. 
dimensionless wavevector. The dotted lines are the type-I so- 
lutions at F = 24 and the solid lines are the type III solutions 
at F = 25. 

The phase below the transition described by the type- 
Ill solutions is characterized by the tendency there of 
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the pancake vortices to stack above each other along the 
field direction. This is clearly demonstrated by the be- 
havior of the pair correlation function hn{r). As can be 
seen from Fig. 3(c) the type-III solutions for /i„(r) for 
n ^ exhibit a singular behavior for small r. This is 
in contrast to the type-I solutions shown in Fig. 3(b). 
The type-I solutions are regular at r = 0. The singular 
behavior at small r in type-III solutions is also reflected 
in the very slow decay at large k of the Fourier trans- 
formed quantities as seen in Fig. 2(b)-(d). We estimate 
that hn{r) ~ a„/r" as r — > where the constants a„ 
and the non-universal exponent a depend on the tem- 
perature parameter F. From Eq. (^ we can also deduce 
that the singularity in /i„ must be cancelled by the same 
singularity in c„, and this is confirmed in our numeri- 
cal calculation. The singular behavior indicates that be- 
low the transition the pancake vortices in different layers 
have a large probability of being on top of each other 
and so behave more like stiff vortex lines rather than in- 
dividual pancake vortices. One way to see the correlation 
perpendicular to the layers is from the n-dependence of 
hn{r) for fixed r. For example for r ~ 0, this quan- 
tity decreases very rapidly to zero as seen from Fig. 3(b) 
for the type-I solutions. This shows that the correlation 
length along the field direction is short. One the other 
hand, Fig. 3(d) shows that, below the transition, /i„(r) 
decays very slowly with n and appears to approach a 
non-zero value suggesting that there might even be long- 
range order present. As will be demonstrated later in 
this section, this is not the case. In fact there seems to 
be no order parameter that describes correlation across 
the layers which becomes non-vanishing below the tran- 
sition. Thus our transition seems to be best described as 
a liquid-liquid transition for which the low-temperature 
phase has much longer correlations along the field direc- 
tion than the high-temperature phase. 

Using the calculated correlation functions, one can also 
determine other thermodynamic quantities like the ex- 
cess internal energy U°^. The dimensionless internal en- 
ergy per particle, u''^ = U°'^/2ttJN is given by 



I c^^r m„(r)/i„(r). 



The specific heat per pancake excluding the ideal gas 
contribution can be calculated by differentiating the in- 
ternal energy with respect to the temperature: Cy = 
N-'^diV^'/dT) = fcB(dM'='^/dr-i). rpj^^ •j^gj,^ Yig. 4 
shows that the specific heat in our system obtained by 
numerically differentiating u^'^. The specific heat is larger 
in the low temperature phase consistent with general 
features of an ordinary phase transition. The liquid-to- 
liquid phase transition observed in Q is quite unusual in 
the sense that the low temperature phase has a smaller 
specific heat. We believe this is also a consequence of 
disorder which in addition to changing the order of the 



transition might also modify the sign of the specific heat 
jump (7|. 
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FIG. 3. The pair correlation function hn{R) for rj = 0.5 
from the type-I (F = 24) and type-III (F = 25) solutions with 
dimensionless R = r/a. (a) ho{R) from the type-I (dotted 
line) and type-III (solid line) solutions; (b) The type-I h„{R) 
for n = 1,2,3 and 4. The overall magnitude of hnir) gets 
smaller for larger n; (c) The type-III hn{R) for n = 1 and 2. 
(d) The type-III h„{ro) as a function of n for very small fixed 
ro ~ 0.016. 
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FIG. 4. The excess internal energy per pancake as a 
function of F. The type-I (circles), type-II (triangles) and 
type-III (squares) solutions are shown. The inset is the plot 
of the specific heat versus temperature F"^. 

For other values of i], the first-order transition point 
Fc can in principle be determined in the way described 
above. This will then lead to a transition line in the 
physical B-T plane on using Eq. (0). Unfortunately, we 
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found it quite difficult using our iterative technique to 
obtain type-Ill solutions starting from type-II solutions 
especially when rj becomes larger, since they have a very 
different form from the type-II solutions especially in the 
small-r region. In the present analysis we were only able 
to identify a couple of transition points which are shown 
in Fig. 5. It is, however, relatively easy to find the spin- 
odal points where the type-I and type-II solutions cease 
to exist. The collection of these points up to very large 
values 77 are shown in Fig. 5. If there were a critical point 
as in an ordinary liquid-gas phase transition, the two 
spinodal lines would have to meet at the critical point. 
Fig. 5 shows that the two spinodal lines almost run par- 
allel to each other in the large j] region. Therefore within 
our model the phase transition seems to continue into 
both the large and small 77 regime, implying that no crit- 
ical points exist. We also note that the slope of the phase 
transition line is opposite to those of the spinodal lines. 
Since the transition line should lie between two spinodal 
lines, we expect that as we go to the large-77 region the 
transition line will change its direction and show a kind 
of re-entrant behavior in the ?7-r space. 



larger systems. The c-axis correlation function shows 
a similar behavior for larger systems as seen in Fig. 6. 
However, from Fig. 6 one can also see that in the limit 
Ni —> 00 the correlation function will decay to zero as n 
becomes large, implying that there is no long-range order 
in the thermodynamic limit. This analysis also tells us 
that in order to study the true thermodynamic behavior 
of the system, free from finite size effects, we will have 
to consider a very large number of layers. This will be a 
very challenging problem numerically. 
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FIG. 5. The spinodal points Ti (inverted triangles) and 
Tii (triangles). The first-order transition points (crosses) for 
T] = 0.45 and 0.5. 



The liquid phase described by the type-Ill solutions is 
characterized by the inter-layer correlation functions be- 
coming singular at small r and decaying very slowly with 
n (the layer separation) at fixed r. For a system con- 
sisting of iV; = 11 layers, Fig. 3(d) shows the correlation 
function almost becomes a non-zero constant for large 
layer separation n. If this continues to be the case for 
larger systems, then it would imply that there is indeed 
long-range order below the transition in the direction per- 
pendicular to the layers. We therefore investigated the 
c-axis correlation function as done in Fig. 3(d) but for 



FIG. 6. The pair correlation function /i„(ro) for very 
small fixed ro — 0.016 as functions of n. The two pa- 
rameters used are F = 17 and rj = 0.5. Systems consist- 
ing of Ni = 11,13, ...,23 layers are presented. Because 
of the periodicity one only has to consider the values of 
n=l,2,...,(iV,-l)/2. 

One aspect which we did not consider in this anal- 
ysis is the possibility of having crystalline solutions of 
the HNC equations. In order to study this one has to 
start from a HNC functional which allows spatial varia- 
tion of the density. Eventually one has to evaluate and 
compare the free energies of the crystalline solution with 
the liquid state solutions. Whether such solutions exist 
and have a lower free energy than the liquid solution is 
left for future work. Other possibilities exist such as the 
crystallization transition (if any) taking place at a lower 
temperature than the liquid-liquid transition. Evidence 
of two transitions has actually been seen in a simula- 
tion of pancakes interacting with screened interactions 
]l^ . On this scenario one would argue that the presence 
of disorder removes the crystalline phase at large fields 
leaving just the liquid-liquid transition. In principle, be- 
cause different symmetries are involved, a liquid-liquid 
transition followed by a liquid-crystal transition on fur- 
ther cooling is a distinct possibility, although intuitively 
one might expect that the temperature interval between 
the liquid-liquid transition and the onset of crystalliza- 
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tion might be small. Maybe only the presence of disorder 
makes the liquid-liquid transition clearly visible as a dis- 
tinct transition. 

In summary we have considered a model for pancake 
vortices in layered superconductors by deriving a pair 
potential for them. Using the HNC approximation we 
have found that there is a first-order transition within 
the vortex liquid phase without any critical end-points. 
No long-range order associated with the transition could 
be identified. However, in the low-temperature phase the 
pancake vortices stack up to form stiff vortex lines. In 
order to study the rich phase structure exhibited by high 
temperature superconductors in a magnetic field it will 
be necessary to include the effect of disorder into the 
present analysis. This will change our first order transi- 
tion into a continuous transition above some critical field 

0. 
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